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We study the entropy production in a macroscopic nonequilibrium system that undergoes an 
order-disorder phase transition. Entropy production is a characteristic feature of nonequilibrium 
dynamics with broken detailed balance. It is found that the entropy production rate per particle 
vanishes in the disordered phase and becomes positive in the ordered phase following critical scaling 
laws. We derive the scaling relations for associated critical exponents. Our study reveals that a 
nonequilibrium ordered state is sustained at the expense of macroscopic time-reversal symmetry 
breaking with an extensive entropy production while a disordered state costs only a subextensive 
entropy production. 
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Detailed balance is the hallmark of the thermal equi¬ 
librium state. A system is said to obey detailed balance if 
the probability current along any microscopic trajectory 
in the phase ^ace is balanced by that along the time- 
reversed one [l[. Consequently, time-reversal symmetry 
is preserved in thermal equilibrium. 


Thermodynamics of nonequilibrium systems, where 
detailed balance and time-reversal symmetry are broken 
with a positive entropy production, has been attracting 
a lot of interests il. Recent studies have been focused 
on microscopic systems with a few degrees of freedom 
where the effect of thermal fluctuations are strong. Un¬ 
der the framework of stochastic thermodynamics, vari¬ 
ous fluctuation theorems are discovered, which provide 
useful insights on the nature of nonequilibrium fluctu¬ 
ations. Theoretical works foster experimental studies of 
microscopic systems such as molecular motors, nano heat 
engines, biomolecules, and so on mi. 


Macroscopic systems pose an intriguing question on 
the level of irreversibility. Consider a many-particle sys¬ 
tem displaying an order-disorder phase transition whose 
microscopic dynamics does not obey detailed balance. 
Does the broken detailed balance result in time-reversal 
symmetry breaking at the macroscopic level? On the one 
hand, one may expect that entropy productions of each 
particle add up to a macroscopic amount irrespective of 
a macroscopic state. On the other hand, if the system 
is in a disordered phase so that all configurations are 
almost equally likely, then irreversibility may not show 
up on a macroscopic level producing only a subextensive 
amount of entropy. A system in an ordered state has a 
lower entropy than in a disordered state. Then, which 
phase produces more total entropy including the system 
entropy and the environmental entropy? These questions 
lead us to the study of the entropy production in a model 
system undergoing nonequilibrium phase transition. 

In this paper, we investigate the emergence of macro¬ 
scopic irreversibility out of microscopic dynamics with 
broken detailed balance. We find that the total entropy 


production changes its character from being subexten¬ 
sive to being extensive as the system undergoes an order- 
disorder phase transition. The entropy production rate 
per particle exhibits critical scaling laws as an order pa¬ 
rameter does in ordinary critical phenomena, and scaling 
relations among critical exponents are derived. Although 
the results are derived in a specific model system, we ar¬ 
gue that the scaling behaviors should be valid for general 
nonequilibrium systems. 

As a nonequilibrium model, we adopt the particle sys¬ 
tem in two dimensions introduced in Ref. 1 ^. This 
model describes a flocking phenomenon of passive parti¬ 
cles. In nature a flock of birds and a school of fish display 
a collective motion [l6| . Such a phenomenon has been 
studied with microscopic models consisting of active self- 
propelled particles moving at a constant speed [13 [3. 
Flocking takes place when particles are subject to an in¬ 
teraction that favors alignment of individual velocities to 
the mean direction. 

The model in Ref. [l^ is composed of passive particles 
in the thermal reservoir instead of active particles. It 
consists of N Brownian particles of mass m in a two- 
dimensional plane of size L x L embedded in a ther¬ 
mal reservoir at constant temperature T. The particle 
density is denoted by p = N/L'^. Let = {xii,Xi2) 
and Vi = ^ = {vii,Vi2) be the position and the ve¬ 
locity of a particle i = I,-- - ,N. We will represent 
a configuration of the whole system with a short-hand 
notation Z = {X,V) with X = {xi, X2, ■ ■ ■ ,X]\[} and 
V = {vi,V2, ■ ■ ■ , vn}. The equations of motion are given 
by 


nil ■ 

+ ( 1 ) 


where 7 is the damping coefficient and ^i{t) = 

{^ii{t),^i2{t)) is the thermal noise satisfying 


(6a/)) =0 

(6a/)6b(^0) = ^^kBTSijSabSit - t') 


( 2 ) 
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with the Boltzmann constant ks, which will be set to 
unity hereafter. The velocity aligning force Fi{V) is 
taken to be 


F,(y) = rt), X (/ X 1)0 = r (/ - (/ • v,)v,), (3) 

where T is the interaction strength, Vi = Vi/\vi\ is the 
unit vector, and 


i=i 


N 


(4) 


The vector / points towards the average direction of the 
particles, and its magnitude A = |/| plays a role of the 
order parameter for the collective motion. Note that the 
force Fi is perpendicular to Vi. It does not work on the 
particle but turns the direction of Vi toward f. The 
interaction is infinite-ranged. A short-ranged version of 
the model was studied in Re f. Il9l |. 

Numerical study in Ref. [15j found that the system 
undergoes a phase transition separating a disordered 
phase (r < Tc) and an ordered phase (T > Tc). Near 
r = Tc, the order parameter scales as (A)s ~ (T — Tc)^ 
and the susceptibility y = N{{A‘^)s — {A)l) scales as 
X ^ |r — Tel”''', where {)s denotes the steady-state en¬ 
semble average. The critical exponents are given by 
(5lv ~ 0.491 and ^jv cs 1.02, where v cs 0.94 is the 
correlation length exponent ^ |r — Td”'^) [l^. These 
exponents are compatible with those of the mean field XY 
model [ 2 ^, [2lj . When the interaction is infinite-ranged, 
the correlation volume is more useful than the corre¬ 
lation length Since the model under consideration is 
embedded in the two-dimensional space, the correlation 
volume is given by S^y = ^ and scales as^y ~ |r — Td”'" 
with V = 2y. 

The velocity-dependent force breaks the detailed bal¬ 
ance and the time-reversal symmetry. We quantify the 
amount of the time-reversal symmetry breaking by the 
entropy production. Suppose that the system evolves 
along a stochastic trajectory Z[t\ = {(X(f), V(t))|0 < 
t < r} for a time interval r. Following stochastic thermo¬ 
dynamics 0, the total entropy production AS'tot['2[T]] 
along a given trajectory Z[t\ is determined by the prob¬ 
ability ratio of Z[t\ against its time-reversed trajectory 


ra.iectc 

Z^[t] = {(X(T-t),-F(T-f))|0<t<r} [8|,[2|i6|. 

In our model, the total entropy production is decom¬ 
posed into three terms as 

A^tot [2] = A5eys [2] - ® + A^v [^], (5) 

where A^syg is the change in the Shannon entropy of 
the system, the second term is the Clausius form for the 
entropy change of the heat bath with Q being the heat 
absorbed by the system, and the last term S'v appears 
only in the presence of a velocity-dependent force 261 



FIG. 1. s versus F for several values of N. Inset shows 
the finite size scaling behaviors of s when F is below, equal 
to, and above Fc = 1.976. The dotted (dashed) line has the 
slope —1 (—1/2). 


and is given by 


= ^f2l\tF,{V{t)) O (6) 

In the steady state, the ensemble average of ASsys 
vanishes. The thermodynamic first law reads as AE = 
Q + W where AE is the change in the total energy 
E = and W = Jo work done 

by the force. Since W = 0, Q = AE and its steady state 
average vanishes. Thus, the entropy production rate per 
particle in the steady state is given by 



We have performed numerical simulations. The equa¬ 
tions of motion in o are integrated numerically by us¬ 
ing the time-discretized {At = O.OI) Heun algorithm [ 2 ^. 
We took m = j = p = 2T = 1 in all simulations. Figure[T] 
shows that s displays a characteristic behavior signaling 
a continuous phase transition. As N increases, s ^ 1/N 
for F < Fc while it converges to a finite value for F > Fc. 
We also measure the susceptibility of the entropy pro¬ 
duction that is defined as 


Xs{T,N,t) 


I 

tN 




(A^v)^ 


( 8 ) 


where AS'v denotes the entropy production of N particles 
in a time interval r. Figured (a) shows the susceptibility 
measured at fixed r = 64. It has a sharp peak at F = Fc, 
which also reminds us of a continuous phase transition. 
The threshold Fc ~ 1.976 is close to the onset of the 
collective motion reported in Ref. [l^ . We will show that 
the entropy production indeed exhibits the continuous 
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FIG. 2. (a) Susceptibility Xs as a function of F. (b) Xs 

versus (F — F^) in the log-log scale. The dashed line has the 
slope —1. (c) Xs versus N for F < Fc. The dashed line has 
the slope —1. 



FIG. 3. (Color online) Scaling plot of versus |F — 

FcIA^^^*^ according to (1141) . The solid (dashed) line has slope 
1 (- 1 ). 


phase transition and that the phase transition is triggered 
by the onset of the collection motion. 

The entropy production can be related to the order 
parameter A for the collective motion. Using the equa¬ 
tions of motion for dvi/dt, the entropy production in ([5]) 
is written as 


27| 


N 


AS^=J2 / dt 


.^0 


i=l 


1 


1 


7 i m 


N 


(9) 




where Vt,. denotes the gradient operator with respec¬ 
tive to Vi and dWi{t) = dt'^i{t'). The last 

term contributes neither to the ensemble average nor to 
the susceptibility because it is of the order of 0(r^/^) 
with zero mean while the others scale linearly with r. 
Hence, it will be ignored. We then introduce the po¬ 
lar coordinate so that the velocity vector is written as 
Vi = {vi cos 9i, Vi sin 9i). The relation dH) for the vector 
f = (Acos^, Asin'f/j) is written as 


1 ^ 
= - y 

N ^ 


i=i 


By using ([3]) and (ITUl) . we can show that [27 1 


N pT 

A5v = y / dt [A, -B, + Qj + 
i=l “'o 


( 10 ) 


( 11 ) 


where Ai = 1^ sin^(V' - 9i), Bi = cos(i/i - 9,), and 
C = r 

^ Nmvi ' 

The expression in (fTTll gives a hint on the scaling be¬ 
havior of the entropy production. The macroscopic vari¬ 
ables A and tp fluctuate much slower than the microscopic 


variables u^’s and 9Ps. Thus, in taking the ensemble- 
average of (fTT]) , we can use the adiabatic approxima¬ 


tion |2l| to replace A^ and A with their ensemble aver¬ 
aged values. Power counting combined with the adia¬ 
batic approximation leads to the conclusion that the en¬ 
tropy production rate per particle scales as s ~ (A^)s 
(A)g (from Ai and Bi) with the 0{N~^) correction (from 
Ci). Therefore, we expect that the entropy production 
rate per particle exhibits a critical power law scaling 


s ^ (T - r. 


( 12 ) 


with the critical exponent 

Pe = 2/3 (13) 

for r > Pc and s ~ 1/7V for P < P^. When N is finite, 
following the standard finite-size-scaling (FSS) ansatz, 
we expect that 

g = jV-/3e/p$ |'(p _ pjjvi/py (14) 


The scaling function $(a;) has the limiting behaviors 
$(a;) ensuring (1131) and $(a;) guar¬ 

anteeing the scaling in the disordered phase. 

The numerical data in Fig. [T] are analyzed according 
to the FSS form with the mean field critical exponents 
/3e = 1 and V = 2. As shown in Fig. [31 the data collapse 
and the limiting behaviors of the scaling function confirm 
the scaling relation in (fT31) and the FSS form of (1131) . 

The total entropy production AS'v is given by the spa¬ 
tial and temporal sum of the fluctuating local entropy 
production rates. We can derive the scaling form for 
the susceptibility Xs in the following way: Near the crit¬ 
ical point, the correlation volume and time diverge as 
~ |r —Fcl”*^ and ~ [h —Fd”*^*, respectively. When 
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FIG. 4. (Color online) (a) ;)j;s(rc, Af, r) versus r at several 
values of N. (b) Scaling plot of Xs(rc, A^, r)A^ versus 
tN-\ 


N and r in the ordered phase (F > Fc), the 

total entropy production AS'v is the sum of the contri¬ 
butions from M = TN/{^t^v) space-time blocks. All the 
blocks are independent because they are beyond the cor¬ 
relation volume and time. Therefore, the susceptibility 
should scale as Xs ^ x M x ^ which 

leads to the scaling form 

X.^(F-Fe)-^« (15) 

with the susceptibility exponent 

7e = t't + - 2/3e- (16) 


This is the hyperscaling relation extended to the systems 
with anisotropic scaling m, 11^. At the critical point, 
the finite-size effect dominates so that 


Xs(Fc, Af,T) 


, T < A® 


(17) 


with z = i/j/F. In the disordered phase, the entropy 
production rate per particle vanishes as s ^ 1/lV, so does 
the susceptibility Xs 1/N. 

The numerical data support the scaling theory. Fig¬ 
ure [2] (b) shows the susceptibility follows the power law 
of (fT51) with 7 e = 1. This exponent value satisfies the 
hyperscaling relation in (I16|) with = 1, 1? = 2, and 
/3e = 1. The 1/N scaling inside the disordered phase is 
also checked in Figure [2] (c). The FSS behavior at the 
critical point F = Fc is examined in Fig. |4l At a given 
N, Xs(rc,lV, t) increases algebraically with r and sat¬ 
urates to a limiting value (see Fig. SKa)). The scaling 
plot in Fig. BKb) confirms the scaling behavior of (1171) for 
T <C and r > N\ 

We have shown that the broken detailed balance leads 
to the macroscopic entropy production only in the or¬ 
dered phase using the analytic scaling theory and the 


numerical simulations. The entropy production per par¬ 
ticle per unit time s is positive but vanishes as 1 /N in the 
disordered phase, while it is hnite in the ordered phase 
following the power law (see (fT^ ). The susceptibility Xs 
vanishes in the disordered phase and follows the power 
law (see (fTSll ) in the ordered phase. The critical expo¬ 
nents satisfy the scaling relations in (IT^ and (fTHl) . 

The quadratic relation s (A)g is crucial in deriving 
the scaling theory. This relation is derived in a model 
system that has a mean field nature. We argue that the 
scaling behaviors are universal in general thermal sys¬ 
tems undergoing a nonequilibrium phase transition be¬ 
tween a disordered phase and an ordered phase. Collec¬ 
tive motions in the ordered phase are characterized by 
the thermodynamic currents Ji of e.g., energy and parti¬ 
cle. The currents are small near the critical point. Thus, 
following the linear irreversible thermodynamics of On- 
sager [3l|, one can assume that Ji = where 

Xj ’s are the thermodynamic forces and Ly’s are the On- 
sager coefficients. The entropy production rate is then 
given hy dS/dt = Y,iXiJi = which 

supports the validity of the quadratic relation between 
the entropy production rate and the current density. In 
stochastic thermodynamics, the total entropy production 
rate is written as the configuration space average of the 
probability current density squared [SI], which also sup¬ 
ports the relation. It would be interesting to investigate 
the scaling relations in m and p6l) in systems with a 
short-ranged interaction. 


The result that the ordered phase costs more envi¬ 
ronmental entropy production may be understood in the 
framework of the thermodynamic second law. Suppose 
that one changes a coupling constant of a system so that 
it relaxes from a disordered phase to an ordered phase 
in a characteristic relaxation time Aeiax- During the pro¬ 
cess, the system entropy decreases at the rate dSsys/dt ^ 
AtSsys/tj-eiax — (tSgys(ordered) iSsyg(disordered))/ti.eiax ^ 
0. The thermodynamic second law requires that the 
entropy production rate should be nonnegative at any 
moment. Therefore, during the relaxation process, the 
environmental entropy production rate should satisfy 
dSejiv/dt ^ dSsys/dt ^ IAtSgygl/treiax; which gives a 
lower bound for the environmental entropy production 
rate. It should be investigated further whether the in¬ 
equality is working in the steady state. We leave it for 
future work. 


This work was supported by the Basic Science 
Research Program through the NRF Grant No. 
2013RIA2A2A05006776. 


[1] C. Gardiner, Stochastic Methods, A Handbook for the 
Natural and Social Sciences (Springer, New York, 2010), 
4th ed. 






5 


[2] D. J. Evans, E. G. D. Cohen, and G. P. Morriss, Phys. 
Rev. Lett. 71, 2401 (1993). 

[3] G. Gallavotti and E. Cohen, Phys. Rev. Lett. 74, 2694 
(1995). 

[4] C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997). 

[5] K. Sekimoto, Prog. Theor. Phys. Supp. 130, 17 (1998). 

[6] J. Lebowitz and H. Spohn, J. Stat. Phys. 95, 333 (1999). 

[7] G. Crooks, Phys. Rev. E 60, 2721 (1999). 

[8] U. Seifert, Phys. Rev. Lett. 95, 040602 (2005). 

[9] G. Wang, E. Sevick, E. Mittag, D. Searles, and D. Evans, 
Phys. Rev. Lett. 89, 050601 (2002). 

[10] D. Carberry, J. Reid, G. Wang, E. Sevick, D. Searles, 
and D. Evans, Phys. Rev. Lett. 92, 140601 (2004). 

[11] D. Collin, F. Ritort, C. Jarzynski, S. Smith, 1. Tinoco, 
and C. Bustamante, Nature 437, 231 (2005). 

[12] V. Blickle and C. Bechinger, Nature Physics 8, 143 

( 2011 ). 

[13] J. Gomez-Solano, A. Petrosyan, and S. Ciliberto, Phys. 
Rev. Lett. 106, 200602 (2011). 

[14] D. Y. Lee, C. Kwon, and H. K. Pak, Phys. Rev. Lett. 
114, 060603 (2015). 

[15] F. J. Sevilla, V. Dossetti, and A. Heiblum-Robles, J. Stat. 
Mech.: Theor. Exp. 2014, P12025 (2014). 

[16] T. Vicsek and A. Zafeiris, Phys. Rep. 517, 71 (2012). 


[17] T. Vicsek, A. Czirok, E. B. Jacob, I. Cohen, and 
O. Shochet, Phys. Rev. Lett. 75, 1226 (1995). 

[18] G. Gregoire and H. Chate, Phys. Rev. Lett. 92, 025702 
(2004). 

[19] V. Dossetti and F. J. Sevilla, arXiv (2014), 1410.3187vl. 

[20] B. J. Kim, H. Hong, P. Holme, G. S. Jeon, P. Minnhagen, 
and M. Y. Choi, Phys. Rev. E 64, 056135 (2001). 

[21] S.-I. Sasa, New J. Phys. 17, 045024 (2015). 

[22] R. Spinney and 1. J. Ford, Phys. Rev. Lett. 108, 170603 

( 2012 ). 

[23] H. K. Lee, C. Kwon, and H. Park, Phys. Rev. Lett. 110, 
050602 (2013). 

[24] C. Ganguly and D. Chaudhuri, Phys. Rev. E 88, 032102 
(2013). 

[25] D. Chaudhuri, Phys. Rev. E 90, 022131 (2014). 

[26] C. Kwon, J. Yeo, H. K. Lee, and H. Park, arXiv (2015), 
1506.02339vl. 

[27] See Appendix for details. 

[28] A. Greiner, W. Strittmatter, and J. Honerkamp, J. Stat. 
Phys. 51, 95 (1988). 

[29] H. Hong, M. Ha, and H. Park, Phys. Rev. Lett. 98, 
258701 (2007). 

[30] M. Henkel and U. Schollwock, J. Phys. A 34, 3333 (2001). 

[31] L. Onsager, Phys. Rev. 37, 405 (1931). 


Appendix: Total entropy production 

It is straightforward to decide whether a deterministic dynamics is reversible or not. Suppose that a system evolves 
from a configuration Z(0) = (A(0),V(0)) to Z{t) = {X{t),V{t)) along a trajectory Z[t] = {Z{t)\Q <t< r}. If 
one flips the velocity in the final configuration and takes the resulting configuration (A(t), —V (r)) as an initial state, 
then the reversible dynamics lets the system follow the time-reversed trajectory Z^[t] = {Z^{t)\0 < t < t} with 
Z^it) = iXiT-t),-ViT-t)). 

Generalizing this idea to stochastic systems, one can define the irreversibility or the entropy production by comparing 
the probability of trajectories Z[t] and Z^[t]. The probability distribution function (PDF) of a given trajectory Z[t] 
is given by P[Z[t]] = 11 [Z[t]; Z(0)] po{Z{0)), where po{Z) is an initial PDF of being in a configuration Z at time 
t = 0 and n[Z[T]; Z{Q)] is a conditional probability distribution of Z[t\ to a given initial configuration Z{Q). The PDF 
for a time-reversed trajectory Z^\t\ is similarly given by P[Z^[t]] = 11 [Z-^jr]; Z^(0)] Pr{Z{T)), where Pr{Z) is the 
PDF at time r which has evolved from no(^)- According to stochastic thermodynamics, the total entropy production 
for a given trajectory Z[t] is given by [8| 


n[Z[r];Z(0)] poiZjO)) 
n[Z«[r];Z«(0)]p.(Z(T)) • 

It consists of two parts as A^tot = A^sys + ASenv, where 

ASsys = - InpriZ^T)) -I- lnpo(2'(0)) 


(18) 


(19) 


is the system entropy change and the remaining term AS'env is the environmental entropy production. 

The environmental entropy production can be written in terms of physical quantities. This task has been done in 
a recent preprint for systems with an arbitrary velocity-dependent force. We make use of Eq. (11) of Ref. [2^ to 
obtain that 


N 


N 


A5'env = -yX!/ 4-— ^ / Fi{V{t))odv^{t), 

i=i -^0 'y Jo 


( 20 ) 


where the notation A{t) odB(t) = . [B{t+dt) — B{t)) stands for the stochastic integral in the Stratonovich 
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sense [i|. Using mdvi = (rndvi — Fidt) + Fidt in the first term, one can further decompose AS'env as 


AS'env = - 


Y.i=i o [mdvi - Fidt) Y.i=i lo ° ^idt 


N 


T 


T 


iT 


Y, / F,{V)odv, 


2=1 ' 


( 21 ) 


The Langevin equation indicates that 

^ /*T / ^ f'^ 

Q = Y ° ~ = ^ j ° dt ( 22 ) 

i=i -^0 V / Jo 

is the work done by the heat bath through the damping force and the random force, namely the heat absorbed by the 
system from the heat bath. The second term is identically zero since Vi -L Fi. The third term is A5'v. This completes 
the derivation of Eqs. (5) and (6) of the main text. In AS'tot, (ASsys — Q/T) is generic in all thermal systems, while 
the others appear only in the presence of velocity-dependent forces. 

Note that the force Fi does not work {W = 0). Consequently, the thermodynamic first law is written as AU = Q, 
where AU is the change in the total kinetic energy E = X)i=i 


Appendix: Derivation of Eqs. (9) and (11) 

The Stratonovich product Fi o dvt is defined as [H 

^ ^ F,{V{t + dt))+F,{V{t)) ^ n-., 

Fi o dvi = - - - dv^{t) = Y Fia{t)dvia + 2^ -g;^dviadvjb + o{dt) , 

a j,a,b 


(23) 


where i,j = 1, • • • , N are particle indices and a, 6 = 1, 2 are Cartesian coordinate indices. We now use the Langevin 


equation to replace mdvi = Fidt — jVidt + dWi, where dW) = dt'^i(t') satisfying that (dWia) = 0 and 

{dWiadWjb) = 2jTSijSabdt. Inserting this into (1^ . we obtain that 


mFi o dvi = \Fi\'^dt + Fi • dWi + -gY^dWiadWjb + o{dt). 

j,a,b 

Since dWia's are independent of each other, one can replace {dWiadWjb) with {2'^T5ij5abdt) [l|. This yields 


(24) 


N 


AS'v = XI / 


JO 


i=l 


— |F,|2 + -V„. -F, 

'jl m 


N 


— X / dW., 


(25) 


which is Eq. (9) of the main text. As explained in the main text, the last term can be neglected. 

The expression for ASv becomes simpler in the polar coordinate. Let Vi and 9i are the magnitude and the polar 
angle of Vi , respectively. The magnitude A and the polar angle ip oi f are given by Ae®^ = ^ . The force 

Fi = T{f — (/ ■ Vi)vi) corresponds to the projection of f in the normal direction of Vi. Thus, one can write 


F = rAsin(i/; - 6i)6i, 


(26) 


where Oi is the unit vector in the polar angle direction of Vi. It is evident that |F| = TAj sin(^/: — di)\- The divergence 
is given by 


N 


V„. - Fi = — ^rAsin(V'-6l,) = — 7^4 X “ ^®) 


Vt dOi 

r I 

Vi N 


j=o 


Xcos(% -6»0 = 


j¥=i 


— Acos{ifj — Oi) 


(27) 


Inserting the magnitude and the divergence of F into (1251) . we obtain Eq. (II) in the main text. 

I- 












